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Plasmonic nanoparticles exhibit optical field enhancement when localized surface plas- 
mon polariton (SPP) resonances are excited [1, 2], The strength of the enhancement de¬ 
pends sensitively on the nanoparticle’s environment and geometry [3]. The enhancement 
is a vital part of phenomena such as surface-enhanced Raman scattering [4], and finds 
application in areas such as biosensing [5], monitoring lipid membranes [6], modifying 
molecular fluorescence [7] and materials characterization [8]. Localized SPP resonances 
occur because of the way the free conduction electrons in metal particles respond to 
light. For many metals at optical frequencies their response is such that the permittiv¬ 
ity is negative - a critical requirement if the nanoparticle is to support a plasmon mode. 

However, metals are not the only materials to exhibit negative permittivity; mate¬ 
rials doped with excitonic organic dye molecules are of interest for photonics [9] and 
may also possess negative permittivity over a small frequency range [10]. Interest 
in such materials as a means to support surface exciton-polariton (SEP) resonances 
has recently been rekindled [11, 12, 13]. An example of this class of material is a 
polymer doped with dye molecules. In a previous work we showed, through experi¬ 
ment and with the aid of a classical model, that polyvinyl alcohol (PVA) doped with 
TDBC molecules {5,6-dichloro-2-[[5,6-dichloro-l-ethyl-3-(4-sulphohutyl)-henzimidazol- 
2-ylidene]-propenyl]-l-ethyl-3-(4-sulphobutyl)-benzimidazolium hydroxide, sodium salt, 
inner salt) may support localized surface exciton-polariton modes. We extracted the 
complex permittivity s{uj) of this material from reflectance and transmittance measure¬ 
ments of thin films using a Fresnel approach [14]. TDBC was chosen because of its 
tendency to form J-aggregates: this leads to a narrowing of the optical resonance, mak¬ 
ing them interesting for strong coupling [15, 16, 17]. More importantly in the present 
context, at sufficiently high concentrations materials doped with such molecules exhibit 
a negative permittivity; it is this negative permittivity that enables these materials 
to support localized resonances. In our previous work [12] a two-oscillator Lorentz 
model [18, 19] was used to calculate the electric field enhancement and field confine¬ 
ment around the nanoparticle supporting the resonance by use of Mie theory [20, 21]. 
The field enhancement and confinement we predicted compared favorably with respect 
to gold nanospheres, albeit over a much narrower spectral range. Here we extend that 
earlier work by going beyond a simple classical Lorentz oscillator model. In doing so, we 
are able to explore new transient phenomena and develop a richer microscopic physical 
picture for the system. 

In what follows we first outline the elements and assumptions of the quantum model 
we have used. We then compute the relative permittivity, £( 0 .;), using this model and 
compare our results with those obtained from a classical model. We next use our model 
to investigate the steady-state response of nanospheres of possessing this relative per¬ 
mittivity, with a focus on the nature of the localized surface exciton-polariton (LSEP) 
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mode. The response of the same particle to a suddenly turned on applied optical held 
is then explored, and the occurrence of transient LSEP modes discussed. 


2. Theory 

The key difference between the work we report here and previous work based on a bulk, 
macroscopic approach [12] is that here we develop an effective medium description from 
a quantum model of the relative permittivity e = 1 + y, where y is the susceptibility. To 
do so we assume that the molecules in our material can be represented as an ensemble of 
two-level quantum systems. For a general material, an applied electric held E induces 
a polarization in the material P oc {d), where (d) is the average dipole moment of each 
molecule (quantum system). Assuming linearity, P is a linear function [22] of E, given 
by: 

P = ^e,E{xe-^-^ + y*e*"*) = N{d) (1) 

where N is the number density of quantum emitters, and E and d are generally time 
and frequency dependent. To hnd y and hence e we need to hnd {d). If we adopt a 
quantum picture, then (d) becomes the expectation value of the dipole moment and 
can be computed from the trace of pd, where p is the density matrix of the system and 
d is the transition dipole matrix. The density operator p is dehned as p = 
where pk are the relative probabilities of hnding a system element in state \k). In order 
to hnd p, a Hamiltonian that describes the system must be determined. In general, the 
Hamiltonian for an open quantum system can be expressed as [23, 24], 

H = Ho+ Hb + Hi, (2) 

where Hq is the Hamiltonian of the isolated system, Hb describes the interaction of Hq 
with the bath, and Hj describes the interaction of Hq with the applied electric held. 

For TDBC molecules in a PVA host medium, Hb should represent the Sum—Q = 129 
intramolecular [25] vibrational modes (where Um is the number of atoms per molecule) 
with a multitude of intermolecular modes. These vibrational modes are responsible for 
induced decay and dephasing in the system [26, 27], along with a small shift in the 
excited state energy of the molecules [28]. Rather than determining Hb directly we 
have made a commonly-used simplihcation, that of incorporating the ehects of the bath 
(vibrationally induced decay and dephasing) phenomenologically by application of the 
dissipative Lindblad superoperator (see below) and by making the assumption that the 
small energy shift can be ignored [23, 24, 29]. 
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For an ensemble of n two-level emitters (molecules), Hq can be written as [30, 


31], 


Hq — /iC(;o|0)(0| 

i=l 


( \ 

n 

fhOi ^|lj)(li| + ) 

j = l 

V / 


( 3 ) 


where |0) is the ground state of the nanoparticle, and |li) represents a single exciton 
excited in the nanoparticle, localized on molecule i, with the other molecules in their 
ground states i.e. |lj) = |0i,..., Ij, In this way, only a single exciton is permitted 

within the ensemble at any time. The hrst term in brackets in Equation (3), 
represents the average energy eigenvalue of a non-interacting molecule in the excited 
state (an exciton). The second term corresponds to inter-molecular coupling, with 
coupling energy Jij. The coupling is taken to be Forster (dipole-dipole) coupling [32, 30] 
since we assume that the overlap between the wave functions of each site are small. The 
corresponding interaction Hamiltonian Hj modeled in the Schrodinger picture [33] and 
written in the same basis is, 

n 

i=l 

where the coupling strength of the dipole to the external optical driving held is dehned 
as Pi = -E{ri) ■ fii, where is the exciton dipole moment. 


Although thorough, a density matrix formed using Equations (3) & (4) would have 
dimension (n-|-l) x (n-|-l), where n is the number of molecules in the system. Given that 
n can be several thousand for even a moderately-doped 100 nm diameter nanosphere, 
solving for such a large matrix would be computationally very demanding, despite con¬ 
sidering only a single exciton in the ensemble. We therefore seek a simpler Hamiltonian 
for a TDBC-doped nanosphere which approximates the formalism above. 


As a hrst step in this process, we identify that for an ensemble of aggregates (where the 
monomers within each aggregate are aligned with each other), the intra-aggregate cou¬ 
pling terms dominate [34]; this enables us to neglect the inter-aggregate coupling terms. 
By making this approximation, our approach to describe a nanoparticle doped with ran¬ 
domly distributed and randomly oriented aggregates is to hrst describe a Hamiltonian 
for a single aggregate, and then to take an orientational average. 


The next step is to note that for a single aggregate, nearest-neighbour couplings 
dominate. The Hamiltonian matrix obtained under this approximation using Equation 
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( huo 0 0 0 ■ ■ ■ 0 \ 

0 J 0 ■ ■ ■ 0 

0 J J ■ ■ ■ 0 

0 0 J ■■■ 0 


V 0 0 0 0 


J 


( 5 ) 


where J is the nearest-neighbor interaction energy. The eigenvalues and eigenstates for 
this Hamiltonian matrix are derived in our Supporting Information. The hrst eigenstate 
is the ground state |0), with energy eigenvalue huo- The second is a set of excited states 
where a single exciton is delocalized over the aggregate [35, 29], 


m) 



n 



jmiT \ 
n + l) 


iij>. 


( 6 ) 


where 1 < m < n. The single exciton transition dipole moment of the aggregate dQi{m) 
for mode m is related to the transition dipole moment of the monomers (assuming 
identical dipole moments) by [36], 


I / N /l — (—1)"^ / nm \ 

This implies that <ioi(m) is zero for even values of m. Even for very modest aggregates, 
~ n > 6, the leading eigenstate (|1)) gives rise to a transition dipole moment a factor of 
three stronger than the next eigenstate, i.e. the leading eigenstate is the ’brightest’ [37]. 
We can take the transition dipole moment of the aggregate as dgi ~ <ioi(l), and use the 
two states |0) and |1) as an approximation for the aggregate, where |1), using equation 
6, is given by, 

|l> = \/^E=“(^)|li>- (8) 

j=i ^ / 


The eigenvalue of |1) is (c./. Supporting Information, Equation S6), 


hzoi = — 2J cos f—-— ) . (9) 

+ 1 / 

This allows us to write hioi = A. The excitation energy of the aggregate is 

shifted from the monomer value by A, and this shift arises from the interaction with 
other molecules in the aggregate. This energy shift has been observed elsewhere for 
aggregates [38, 28], and is loosely termed the ‘effect of aggregation’ [32]. The magnitude 
of A is typically hundreds of meV [30]. Therefore, by considering only the ground state 
and this (brightest) excited state. Equation (3) can be re-written as. 
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Hq ~ /iC(;o|0)(0| + {hw[ ^ + A)|l)(l| 

= nwo\o){o\ + hooi\i){i\. (10) 

The interaction Hamiltonian for the aggregate is written as, 

^z = G(|0)(l| + |l)(0|), (11) 

where G = E -d^i is the coupling strength of the electric held, -E, to the dipole moment, 
doi, of the aggregate. The Hamiltonian formed by adding Equation (10) & (11) can 
be applied to an ensemble of randomly-distributed aggregates by taking G = E ■ d^i, 
where doi = doi/V is the orientational average of the aggregate dipole moments in the 

system of interest. Here, V is equal to either 2 or 3, corresponding to the number of 

spatial dimensions in the planar and bulk cases respectively (this orientational average 
is derived in our Supporting Information: see section 3). 

Our goal now is to hnd an effective medium value of e at time t and at the frequency 
of illumination, u. The hrst step is to note that Hj/\E{t)\ dehnes the transition dipole 
matrix d for the system as a whole. Given that the density matrix (p) can be used 
to obtain the expectation value of an observable, we seek p using the Liouville-von 
Neumann equation [39], 

p{t,u) = - ^-[H,p{t,uj)] +LDp{t,uj). (12) 

The hrst term in Equation (12) governs unitary evolution. The Lindblad dissipation 
superoperator [40, 41] Ld-, is used to account for the decay and dephasing ehects the 
bath has on the system. In this work, we assume the electron-phonon coupling to 
be weak at room temperature and for weak helds, and this enables the Born-Markov 
approximation upon which this formalism relies [42] to be used. The total dephasing 
rate of the transition |0) <H- |1) is Toi. This quantity is related to the population decay 
rate for the |1) —)■ |0) decay channel, 701 , and the pure dephasing rate, Tg^^^ by [40], 

r„, = ^ + 77 (13) 

The pure dephasing rate, Tg'i^^ arises from phase-changing interactions of the excitons 
with the environment [43], i.e. the bath of vibrational modes. A less approximate ap¬ 
proach could be adopted [42, 44], but assuming a simple rate for Tg']]^ is sufficient for our 
present purposes, that of enabling an illustrative calculation to be carried out. 

For our two-level system. Equation (12) is used to hnd 4 coupled differential 
equations [45], the well-known Optical Bloch Equations [46] (OBEs). Solving these 
for the applied cosine potential allows us to use the Rotating Wave Approximation 
(RWA) [47], as detailed in our Supporting Information. To derive an expression for the 
permittivity, {d) is determined using (d) = Tr(pd). For our aggregates, this value is 
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equal to poi<^oi + c.c.. By choosing the forward-propagating electric held, Equation (1) 
is re-arranged to give e for an ensemble of two-level molecules as, 

e(f, u) = eb+— (14) 

^0 1^1 

This expression is applicable to an ensemble of molecules with number density N, 
arranged in aggregates, distributed randomly spatially and orientationally (in two or 
three dimensions), in a medium of background permittivity e = £& in the single-exciton 
regime. This formula holds for weak helds, as we show below. At hrst glance. Equation 
(14) may appear to diverge for the case where \E\ —)■ 0. The resolution to this is that 
Poi is linear in \E\ in this limit; this gives us a held independent permittivity in this 
limit as expected. 

3. Results and Discussion 

Eor our model we require the following parameters: doi, hcui, 701 and Tq)!^. We used our 
experimental rehectivity and transmittance data (for a 1.46 wt% TDBC:PVA 70 nm 
him [12]) to determine that hui = 2.11eE = 588 nm. This agrees with the values 
obtained by van Burgel [48] and Valleau [30] although it is a slight change from our pre¬ 
vious work, where we indicated that the transition occurred at 2.10 eV (590 nm), with 
a (weaker) shoulder transition at 2.03 eV (610 nm). Our revised value follows from an 
improved Kramers-Kronig analysis of our original data, as outlined in our Supporting 
Information. 

From photoluminescence measurements [49], we took the decay rate of |1) to be 
7oi = 1.15 X 10^^s“^ for the aggregate in a PVA host medium. Using Molinspiration ©, 
we determined the molecular weight and the ehective volume of the TDBC molecule. 
Together with the concentration of the solution, these quantities allowed us to determine 
the molecular number density to be A^ = 1.47 x 10^^ m~^. We were then able to esti¬ 
mate the transition dipole moment for TDBC molecules in aggregate form, {d^i), and 
the dephasing rate, Pq([\ by htting the steady-state solutions for Equation (14) to our 
experimental data for e(a;) by adjusting |<ioi| and Pq([^ In this way we found the dipole 
moment to be 48 debye (D). The TDBC-doped thin Elms from which the experimental 
data were obtained were produced by spin-coating [12]. Previous work to investigate 
the orientation of dipole moments in thin polymer films produced by spin-coating found 
that the dipole moments he predominantly in the plane [50]. Assuming that the TDBC 
aggregates also lie in the plane of the spun films reported in [12], then the value of 48 
D we have determined here is a two-dimensionally averaged value, implying that the 
on-axis dipole moment of an aggregate is doi = 97 D, and the three-dimensionally av¬ 
eraged moment is 32 D. This three-dimensionally averaged moment compares with the 
24 D estimated by van Burgel et al. [48] from experiments in solution (3-dimensional). 
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The dephasing rate, Tq'^^ was found to be equal to 17 meV, which is ~ kbT as 
expected [30]. To provide additional support for our value of Tq'^^ we extracted 
and modeled e{u) from the reflectance and transmittance data from a 5.1 nm thick 
him obtained by Bradley et al. [51]. We determined Tq*^^ to be around 13 meV, a 
value comparable with our own, bearing in mind that different bath spectral densities 
associated with differences in the host and substrate may change the value of Tq'^^ Our 
results for e{u}) against experimentally-determined data for our him are displayed in 
Figure 1. 



Figure 1. Extracted relative permittivity e(w) from a 1.46 wt% TDBC:PVA film 
(circles for real part and crosses for imaginary part) [12]. The theoretical fit for e{uj) 
indicated by the solid and dashed lines corresponds to a concentration of 1.46 wt% 
(3.22 wt%) for a two (three) dimensional distribution of dipoles. 

The permittivity of the thin him shown in Figure 1 has been modeled assuming the 
dipole moments of the aggregates are randomly oriented in the plane of the him. In 
what follows we wish to look at the optical response of a nanoparticle. For generality, 
and to ensure we consider an isotropic system, we will consider the dipole moments of 
the aggregates to be randomly oriented in three dimensions. Making this assumption 
requires us to increase the number density of our molecules from 1.47 x 10^^ m~^ to 
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3.29 X 10^^ m~^. In this way, our nanoparticle will be comprised of a material that 
has the same permittivity as that shown in Figure 1. In all the calculations that follow 
we use this number density, which corresponds to a concentration of TDBC in PVA of 
3.22 wt%. 


3.1. Numerical Results: Steady-State 

We now explore theoretically the Mie [20, 21] absorption efficiency spectra Qabs{(.^) for 
a 100 nm diameter nanosphere of 3.22% TDBC:PVA, assuming a volume distribution 
of dipole orientations, based on e{u) calculated using Equation (14). In practice, the 
applied optical held we model here might be a laser beam. For a 1 mW laser with a 
spot diameter of 1.5 mm, the strength of the electric held of our incident optical held 
would be equal to 462 Vm~^] we assume this value here. 

In a 100 nm diameter nanosphere of our material, there are on average n = 1.72x 10^ 
molecules. Note that it is the number of molecules and by extension their number den¬ 
sity, which is the important quantity (and is used for N in Equation (14)) rather than 
the number density of aggregates, since each molecule provides a potential site for exci- 
ton excitation. To check the validity of our assumption that multi-exciton and nonlinear 
ehects [43] can be neglected, we computed the maximum expectation value of the num¬ 
ber of excitons in the nanosphere {uex = max(pii)n) using Equations (10) and (11) in 
Equation (12). We found that n^xln 1 holds for laser powers of up to 10^ W with 
a spot size of 1.5 mm. Given that our laser power is 1 mW, we assumed that the 
single-exciton linear regime is sufficient to describe the system under this illumination 
power. 

In Figure 2 we plot the absorption efficiency Qabsi^^) for a 100 nm diameter nanosphere, 
calculated for a variety of permittivities; in each case the absorption efficiency is calcu¬ 
lated using Mie theory [20, 21]. Calculated values for Qabs based upon the permittivity 
obtained using our improved analysis of experimental data are shown in Figure 2 as 
a dashed line. Our quantum theoretical spectrum for Qabsi using e(Lo) from Equation 
(14), is shown as the solid line. This theoretically derived spectrum provides a close 
match to the extracted data, most importantly for energies in the region of interest 
below 2.22 eV. For energies exceeding 2.2 eV, there is a limb in the extracted data 
(dashed curve) which might perhaps be attributed to inhomogeneous (non-Lorentzian) 
broadening which is not accounted for using the OBEs. Also displayed in Figure 2 is 
the result for Qabs using a best-£t classical Lorentz oscillator model (the parameters for 
which can be found in our Supporting Information) shown as a dotted line. It can be 
seen that the quantum model outlined in the present paper provides an improved £t 
to the experimental data. We attribute this to the inclusion of dephasing (Fq*^^) in the 
model: if Fq^ were set to zero, a Lorentz model would be recovered in the steady state. 
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Figure 2. Mie calculations for the absorption efficiency Qabsioi) in the steady-state 
for a 100 nm diameter nanosphere of 3.22 wt% TDBC:PVA using the values for e{uj) 
from experiment (dashed line), from two-level OBEs (solid line), and from our previous 
Lorentz oscillator model (dotted line). The material absorption coefficient k (imaginary 
part of the refractive index), normalized to unity, is also plotted for illustrative purposes 
(long dashed line). Inset: a 3D representation of the aggregated emitters (assuming 
brick-stone aggregation, with 15 molecules per aggregate) randomly distributed in a 
100 nm diameter nanosphere. 


and the single damping term in the Lorentz model would have to accommodate both 
decay and dephasing. Therefore, by including dephasing, the actual physical value of 
the decay rate 701 can be included to achieve an accurate result for £. 

It is interesting to note a key feature shown by the data in Figure 2: Qabs-, reaches its 
peak value at 2.16 eV (574 nm). This is in contrast to the absorption coefficient, k{u:), 
which peaks at 2.12 eV (586 nm), shown as a long dashed line in Figure 2. This differ¬ 
ence in spectral position arises because the peak in Qabs is not due simply to absorption: 
rather, it is due to the excitation of a localized SEP mode [12]. Conhrmation of this 
interpretation comes from two sources. First, in the quasistatic limit the polarizability 
of the nanosphere follows the Clausius-Mossotti condition, for which resonance occurs 
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Figure 3. The time-averaged electric field strength normalized to the incident field 
strength (color plot) and the Poynting vector S (arrows) in the vicinity and on the 
surface of the 100 nm 3.22 wt% TDBC:PVA nanosphere, with incident power flow 
along the positive z-direction. Data are calculated for incident photon energies of a). 
2.16 eV = 574 nm and b). 2.12 eV = 586 nm corresponding to peak absorption 
efficiency and k respectively. 


when e is real-valued and equal to —2 (when the nanosphere is in free space) [52]. From 
Figure 1 this can be seen to be approximately true for our absorbing 100 nm diameter 
nanosphere, as the permittivity value at the wavelength of peak absorption efficiency, 
Qabs, is complex and equal to e = —2.251 -|- 1.728h This difference from e = —2 origi¬ 
nates from the fact that e = —2 only gives the resonance condition if the imaginary part 
of e is zero; the complex nature of the permittivity changes the spectral location of the 
absorption peak. Second, Qabs near the peak goes well above unity: this is associated 
with held enhancement [2], another signature of a resonant mode. The enhanced electric 
held in the vicinity of the nanosphere is illustrated graphically in Figure 3(a), together 
with direction of power how shown by the Poynting vector S. 

In Figure 3, the incident electric held is polarized in the x-direction. The Poynting vector 
arrows shown in the hgure were calculated at starting points for which 2 ; = —200 nm 
and a; = 0 nm, linearly spaced in the range —200 nm < y < 200 nm. Subsequent points 
for evaluation of the Poynting vector were taken at 10 nm steps in the direction of the 
Poynting vector at each point, resulting in the hux lines shown. The power how in 
Figure 3(a) shows that incident light is drawn towards and absorbed by the nanosphere 
for starting positions up to around 130 nm from the central position of the nanosphere. 
This demonstrates that at this energy, the nanosphere absorbs more light than the light 
geometrically striking it [53], and hence Qabs >1- In comparison, absorption at the 
transition energy, i.e. at 2.11 eV is seen only as a shoulder mode in the absorption 
efficiency of the nanosphere (Figure 2) and the efficiency does not exceed unity. The 
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power flow around the nanosphere for the energy at which k peaks (2.12 eV) is shown 
in Figure 3(b), and the enhancement of the held is much weaker than for excitation on 
resonance at 2.16 eV. 

3.2. Numerical Results: Time Domain 

We now turn our attention to the time domain. Our theoretical model for dynamic 
processes in two-level quantum systems subject to a perturbing cosine electric held is 
similar to models considered elsewhere [54, 55, 56, 57], but here the observable of interest 
arises from the temporal evolution of the coherences of the density matrix, rather than 
the populations. The dynamics of a two-level ensemble subject to a pulse potential has 
been the subject of recent investigation [58], but here we investigate a rather diherent 
case: that of a cosine potential of hxed amplitude that is switched instantaneously on 
at some moment in time. We do this to provide an easily soluble model that illustrates 
the time-dependent phenomena we wish to discuss. 



Figure 4. Calculated Qabs{t,uj) for a 100 nm diameter nanosphere of 3.22% 
TDBC:PVA warming up a pure state at t = 0 with a 1 mW laser set at five different 
detunings S from the exciton transition. 
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By using Equation (14) to calculate e(t) for a given illumination frequency u as 
before, Qabs(t,u) can be determined and its temporal behavior examined. To do this, 
we again use Mie theory. This is an approximation since the helds scattered in Mie the¬ 
ory are assumed to be instantaneous. Given that the dynamics seen in Figure 4 evolve 
over a few femtoseconds and that light propagates over a length scale three times the 
size of the nanoparticle during a single femtosecond, and that the nanoparticle is illumi¬ 
nated with an electric held of constant amplitude, this approximation is deemed to hold. 
Mie theory can therefore be used to give a quasi-instantaneous picture of the absorption. 

Qabsit) is shown in Figure 4 for hve different detunings, 6 = fvjji — hca, from the 
transition at tvjJi = 2.11 eV. We assume that all the molecules in the nanoparticle are 
initially in their ground state. At t = 0 we turn on our held abruptly. We see that 
a steady-state response is attained after > 200 /s, but interestingly, for 5 = 0.09 eF, 
Qabsit) repeatedly exceeds unity in spite of the steady-state value of Qabs being below 
unity at this detuning. Since Qabsit) > 1 implies held enhancement, these data are 
indicative of a transient LSEP mode being present at early times. The time-dependent 
behavior comprises two contributions: the hrst is the oscillatory behavior arising from 
Rabi oscillations; the second is the transient ehects associated with the sudden turning- 
on of the held. In the latter, the magnitude of the density matrix coherences exceed 
their steady-state values for tens of femtoseconds, resulting in larger values of e and 
hence diherent absorption properties to the steady-state. 

If e passes through the Clausius-Mossotti condition for the nanoparticle as it 
approaches its steady-state value, Qabs exceeds unity, implying a transient LSEP mode. 
This is seen best for a detuning of 0.091 eV between 0 — 30 fs. The Rabi oscillations 
follow the generalized Rabi frequency Clu, given by [59], 

nR = (15) 

where, Gi 6, and CIr —)■ 5 in this case, Gi and 6 are the Rabi frequency and the de¬ 
tuning respectively. These Rabi oscillations are naturally convoluted with the transient 
ehects. This implies, together with the short timescales involved in the system, that it 
would be a challenge to see these transitory ehects, but might perhaps be possible [60]. 
Critical to the transient LSEP lifetime is Pg^^ If this dephasing could be reduced with¬ 
out losing the transient negative permittivity that is essential for held enhancement (and 
held conhnement), then the transient timescale of the system would be increased up to 
a maximum of l/7oi- This corresponds to the picosecond regime for our TDBC:PVA 
system. Under this circumstance, transient LSEP modes would become more easily 
observable. 
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We have re-evaluated the measurements reported in our previous work and have 
obtained an improved permittivity for our J-aggregate-doped 1.46 wt% TDBC:PVA 
polymer him. Using a quantum-mechanical framework we have given support to our 
previous investigation based on a classical analysis [12], that TDBC doped nanoparticles 
can exhibit a localized surface exciton-polariton (LSEP) mode. We have used a 
quantum model to show that these nanoparticles may also exhibit transient LSEP 
modes in the sub-picosecond regime. These results help strengthen the idea that 
molecular excitonic materials provide an interesting alternative upon which to base 
nanophotonics [9]. By using molecular materials the possibility of bottom-up approaches 
such as supramolecular chemistry and self-assembly can be brought to bear on the 
production of nanophotonic structures. 
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1. Derivation of the excited states of an aggregate 


The eigenvalues of the matrix Hamiltonian in our main text (Equation 5) are solutions 
to the equation, 
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0, (SI) 


where J represents the nearest-neighbour coupling. Equation (SI) can be expressed in 
the form, 


hojQ — A 0 
0 Hn 


(S2) 


where Hn is the sub-matrix written explicitly in Equation (SI), and is representative 
of an aggregate with n molecular units. The eigenvalues of Equation (SI) can be 
determined by solving. 


{nuj^-\)\Hn\=Q. (S3) 

The hrst of these is the ground state energy Aq = huo, which readily yields the 
eigenvector |0). The other eigenvalues and eigenvectors require more consideration. 
The hrst step is to hud an expression for \Hn\- Writing this determinant in terms of 
further sub-matrices gives the following recursive relationship, 

\Hn\ = (M'^ - ^m)\Hn-l\ - j'|Rn-2 


(S4) 
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The next step is to identify that this recursive relationship can be put into the same 
form as the recurrence relation for Chebyshev polynomials [1] i.e. 

Unix) = 2xUn-lix) - Un-2ix). (S5) 

In doing this, and after a little algebra, an expression for the eigenvalues is 
determined, 

Xm = - 2 J cos ( ^ , (S6) 

where 1 < m < n. The hrst of these (Ai) is the excited state of the aggregate (hui) 
taken in our main text. 


The eigenvectors \m) are determined by analysis of The element of 

Hn\m) is. 


Jnij^i + hoo^^mj + Jrrij+i = Xr 


,171 


3^ 


(S7) 


which can also be put into the form of Equation (S5). Making this identihcation and 
a little algebra yields the following expression for the normalised excited states of the 
aggregate [2, 3], 
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i=i 


sm 


( jmn 
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(S8) 


2. Solving the Optical Bloch Equations 

The Optical Bloch Equations (OBEs) derived from the Liouville von-Neumann equation 
can be written in the compact form as, 

p = Lp, (S9) 

where p is a vector of the density matrix elements. One may solve the OBEs by 
application of the Rotating Wave Approximation (RWA) [4] with subsequent use of a 
matrix inversion method, or numerically by way of a Runge-Kutta method, such as the 
RK10(8) method [5]. In this section, both of these approaches are detailed, evaluated, 
and shown to be equivalent for an ensemble illuminated with a cosine potential. 


2.1. Matrix inversion method 

The matrix-inversion method relies on application of the rotating wave approximation 
(RWA) [6] to the perturbing potential, which enables one to write Lit,u) as a time- 
independent matrix, Liu). In this case, solutions for the density matrix elements can 
be written in the following form, 

Pmnit,U3) = (SIO) 
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using the principle of superposition. The coefficients Cj are determined from initial con¬ 
ditions and from the eigenvectors of L{u). The angular frequencies cn* are related to 
the eigenvalues A* of the matrix L{u) by iui = Aj. Aj is complex with a negative real 
part in order to conserve probability. The RWA may be applied to a cosine potential, 
G = d-Eo cos{ut). An advantage of this method is that computation time is very short 
and weak helds may be considered easily. 

A restriction on the RWA is that only a held of one frequency may impinge upon the 
system in this method, since the time-independence of L must be preserved. In order 
to solve the OBEs for circumstances including a frequency-spread pulse, a numerical 
approach must be used instead, as we now indicate. 


2.2. Explicit Runge-Kutta methods 


For the general case where L cannot be written as time-independent e.g. when the 
system is subjected to a pulse, a numerical method must be used to solve the OBEs. 
The OBEs have the general form y = f{t,y) which has a general solution written in 
discretised form. 


Vn+l Vn T hn ^ ^ biki, 


2=1 


where, 

/ i=*-i 

hi y I tn -1“ hnOj, Pn “I” hn ^ ^ ^ijhj 

V j=i 

The estimated error at each step is evaluated as. 


(Sll) 


(S12) 


Cn+l OC y-n+l Un+l- 


(S13) 


The values a*, R and Cij all depend upon the specihc method involved. The original 
Ist-order algorithm using this general method is the Euler method [7], but a much more 
accurate method is the Runge-Kutta (RK4) method, which has been used previously 
to probe the dynamics of 2-level systems [8]. This method suffers from being non- 
adaptive, in the sense that the step size is always taken to be constant. This makes 
numerical solutions for rapidly-changing behaviour unreliable. An improved approach 
is to dynamically allocate the step size between each iteration through a comparison of 
the estimated local error between the 4th and 5th-order solutions at each point. The 
yields the Runge-Kutta-Fehlberg (RK4(5)) method [9]. The advantage of an adaptive 
numerical method such as the RK4(5) method over a non-adaptive method such as the 
RK4 method is that rapidly-changing behaviour can be modelled with greater accu¬ 
racy. Local errors are also minimized and the resultant solution that one determines 
is numerically smoother [10]. However, for our work the RK4(5) method is insufficient 
for producing numerically stable solutions for weak helds. The RK10(8) method (a 
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Figure SI. a). The relative occupancy of the excited state pn, for a 1 mW laser 
potential with spot size 1.5 mm. b). The relative occupancy of the excited state pn, 
together with the coherence poi (real and imaginary parts shown) for a 10 MW laser 
with the same spot size. 


lOth-order Runge-Kutta method with in-built error estimation by comparison to the 
8th order) [5] was tested and found to produce smoother output for insignihcantly more 
computing time than that of the RK4(5) method. 

The main advantage of a Runge-Kutta method over the RWA is that no terms in 
the Hamiltonian are neglected. As a result, any arbitrary potential can be modelled, 
not just a cosine potential. However, the main drawback of any Runge-Kutta method 
is that the computation time for the process is many times longer than using the RWA 
and matrix inversion. This longer timescale is sometimes prohibitive depending on the 
parameters involved. A further drawback is that for weak helds the solutions to the 
OBEs become stiff in time and numerical instability is encountered. This can be seen 
from the fact that with our calculations for TDBC illuminated with a 1 mW cosine 
potential with spot size 1.5 mm, the relative occupancy of the excited state pn does 
not exceed one part in 10^° and so poo remains more-or-less constant in time with value 
1 as shown in Figure SI (conversely, the coherences oscillate in time at the applied 
frequency and are of order 10“®). Therefore, the populations are incredibly stiff and 
this is where explicit numerical methods fail. One way to avoid the resultant numerical 
instability is to conhne this procedure to solving OBEs with strong helds. To achieve 
physically meaningful solutions in this circumstance, one would need to introduce other 
effects such as multi-exciton recombination and nonlinear effects. By neglecting these 
effects, similar solutions can be obtained to those of weak helds for the coherences, pro¬ 
vided the held does not induce signihcant population inversion. As an alternative to 
this compromise, one can use an implicit Runge-Kutta method to solve stih equations. 
However, the computation time for implicit Runge-Kutta methods was found to exceed 
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that of the explicit methods by at least a factor of ten, making this approach ungainly. 

In summary, our procedure is as follows: for weak helds the RWA is used. Solutions 
derived using this method are supported by the solutions obtained using the RK10(8) 
method for strong helds with multi-exciton and nonlinear effects neglected. We conclude 
that by careful use of the RWA, realistic behaviour can be simulated subject to the 
condition that only one frequency illuminates the system. For slowly time-varying 
strong helds or for two strong lasers, an explicit Runge-Kutta method may be used, 
but for weak time-varying helds an implicit Runge-Kutta method may be implemented 
successfully. 

3. The permittivity of a collection of randomly oriented non-interacting 
dipoles 

The quantum mechanical model discussed in the main text predicts the polarizability 
of an individual aggregate. A typical macroscopic sample consists of a large number of 
randomly oriented aggregates that can be treated as non-interacting. The permittivity 
of such a macroscopic sample has a simple relationship to the microscopic aggregate 
polarizability that depends on the dimensionality of the sample, which we derive in this 
section. 

For a collection of N dipoles with dipole moments d* the polarization per unit 
volume P is given by, 

1 N 

i=l 

In the second step of the above equation we assumed that the dipole moments are 
induced by an electric held that can be treated as uniform over the volume V, and 
that the dipoles have a tensor polarizability cXi{u) that depends on the frequency of the 
held. The average value of the polarizability is dehned as (Q:(to')) = A^~^ CKi(ca). If 
we assume that the dipoles are of the same type then the polarizabilities cx.i{u) have the 
same magnitude a{u), and diher only due to the dipole orientation fii. For a collection 
of such dipoles distributed in a P dimensional space. 


(Q:(a;)) = a{uj){fi ® h) 


= a{u) {{nl)x (g) * -h {ny)y <^y 


a{u) 

V 


{x<S)x + y<S)y + ...) 


(S15) 


where to obtain the second and third lines we used isotropy to determine that quantities 
such as (n^Uy) are zero and that (n^) = (n^) = ..., and used the normalization of the 
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unit vectors rij ■ = 1 to determine that (n^) + (n^) + ... = = 1. Combining 

Equations (S14) and (S15) we can then End the macroscopic permittivity £ which is 
isotropic and equal to, 


e{u) 


1 + 


Na{oj) 

V 


(S16) 


with V = 2 for a planar sample, and V = 3 for a bulk one. 


4. Improved Analysis of Experimental Data 

To extract the complex relative permittivity £ or equivalently, the complex refractive 
index n = n + m, of our thin film, we compared our experimental values of 
reflectance {Re) and transmittance (Te) at normal incidence with theoretical ones for 
each wavelength, following the theoretical framework outlined by Heavens [11]. In this 
formalism, the function, 

/(n, k) = \Tt{n, k) - Te\ + \Rt{n, k) - Re\, (S17) 

must be minimised, where the subscript t denotes theoretical values. This process relies 
upon perfect values of Re and Te and n and k are ‘guessed’ to he in a sensible range. 
Eq. S17 is minimised over this range. By this process, small errors in R^, and Tg coupled 
with rounding errors in the guessed range of values for n and k can yield spurious final 
values of n and k. In order to improve this process, we made additional calculations for 
K, varying the thickness of the film from 63 — 77nm, the range of experimental uncer¬ 
tainty. Our results are shown collectively in Figure S2(a). 



Figure S2. k for a lA6wt% TDBC:PVA film, assuming a range of thicknesses before 
(a) and after (b) adjustment for these thicknesses. 

In this figure, we have kept the first two minimum values of Eq. S17 for each wavelength. 
Therefore, the figure shows the physical and first spurious solutions. The gaps associated 
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with taking a single thickness are evidence for the need to consider other thicknesses. 
Doing this produces a range of values for n for each wavelength. We then adjusted each 
value for the thickness taken by assuming for a given thickness t, 




= e 


-K't' 


(S18) 


which leads 


to, 


K = —K. 

t' 


(S19) 


This has the effect of deconvolution on our data, to produce Figure S2(b). The two 
zero-values show where the process has still failed. Our final step was to eliminate the 
spurious values for k, and use the Kramers-Kronig relations to find n, as shown in Fig¬ 
ure S3. 



Figure S3. The improved extracted values of the real (blue circle) and imaginary 
(black pluses) parts of the refractive index from experimental values. 

This process produces much smoother values of h than in our previous work. This is 
because we assumed previously that the true value of n lay in the middle of the two 
solutions found in Figure S2 and the input thickness was held constant. With this im¬ 
proved analysis by varying the input thickness, it can be seen that intermediate values 
between the two clear solutions in Figure S2 do not exist, and it becomes apparent that 
the upper solution is a spurious one. 
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Figure S4. The real (blue) and imaginary (black) parts of the a), permittivity and 
b). complex refractive index of our 1.46 wt% TDBCiPVA film. Solid lines correspond 
to results from Optical Bloch Equations and dashed lines correspond to results from a 
single-oscillator Lorentz model. The discrete data points correspond to our extracted 
values. 


A single Lorentz oscillator model can be nsed to describe the permittivity of a 
material classically [12, 13], 


s{u) — Em + 


Uq — oj'^ — iu'fo 


(S20) 


The parameters for this model which £t our data in Fig. S3 best (assuming a host 
medium of = 1.52) are Uq = 2.11 eV, /o = 0.3, 70 = 46.1 meV. The output of these 
for e and the complex refractive index n are illustrated in Figure S4 together with our 
quantum £t and extracted data. It can be seen that the Lorentz model is not as close 
to the extracted data as the quantum model, particularly for the real parts of both e 
and h. 
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